Benvenute e benvenuti al laboratorio di cosmologia!
Nelle celle di questo notebook esploreremo alcune funzionalità di base ed applicazioni di alcun pacchetti di python sviluppati dai cosmologi per le proprie ricerche. I pacchetti in questione sono a tutti gli effetti delle "calcolatrici scientifiche": così come la nostra calcolatrice può estrarre radici e risolvere equazioni algebriche di base premendo il pulsante uguale, noi, al prezzo di fornire alcune istruzioni in più (dare istruzioni al computer è il significato ultimo di programmare) potremo sfruttare le potenzialità di Python per giocare con (o all') universo.
In particolare, esploreremo le seguenti funzionalità:
Per contatti: andrearubiola.97@gmail.com
La "calcolatrice scientifica" va inizializzata. Importeremo quindi per prima cosa dei pacchetti di utilità generale che consentono di fare i calcoli fondamentali per ogni fisico, nonché eponimi dei primi esami che sempre si danno all'università: algebra (con $\texttt{numpy}$) e analisi $\texttt{scipy}$, che consente anche applicazioni di statistica, interpolazione e ottimizzazione (ricerca di massimi e minimi) di funzioni.
In aggiunta, $\texttt{Matplotlib}$ consente di plottare grafici ed in generale gestire la grafica dei nostri risultati.
######## Utilità generali #############################
import numpy as np # algebra e manipolazione di vettori (array) di dati
import scipy # analisi e manipolazione di funzioni (integrali, derivate, equazioni, interpolazioni...)
######### Grafica #####################################
%matplotlib inline
import matplotlib.pyplot as plt # per la grafica
# Latex: per impostare la grafica "come se fosse un vero articolo scientifico"
font = {'size' : 25, 'family':'STIXGeneral'}
axislabelfontsize='x-large'
plt.rc('font', **font)
plt.rcParams['text.usetex'] = True
np.set_printoptions(precision=3,linewidth=200,suppress=False)
La ricetta generale per fare cosmologia si può ridurre ai seguenti ingredienti:
Le misure sulla struttura a grandi scale (large scale structure, LSS) o sulla radiazione cosmica di fondo (cosmic microwave background, CMB) sono in gran parte statistiche: queste quantità ci consentono di rappresentare, interpretare e fittare i dati per estrarne i parametri che intendiamo misurare. Esistono d'altra parte anche osservabili non statistiche funzione di quantità fisiche di interesse, gran parte delle quali rientrano nelle cosiddette "misure dirette": in particolare, è possibile misurare la costante di Hubble in funzione del redshift, ad partire dalla luminosità intrinseca delle supernove Ia, o, tramite misure di lente gravitazionale forte, "pesare" gli ammassi di galassie e stabilire il valore del parametro $\Omega_{materia}$. Più misure indipendenti abbiamo a disposizione, meglio è, perché aumenta il nostro potere vincolante (constraining power) sul modello, che ci consente di valutarne accuratezza (vicinanza al valore vero di una misura) e precisione (piccole dimensioni delle incertezze).
Un attimo, miei venticinque lettori!
Ora, le quantità statistiche possono essere calcolate nello spazio reale (delle posizioni, delle velocità...), cioè direttamente dai dati osservati. Ad esempio, la funzione di correlazione a due punti è la corrispondente ("trasformata") nello spazio reale dello spettro di potenza ed è un concetto abbastanza intuitivo: è l'eccesso di probabilità di trovare coppie di oggetti -- in particolare, galassie -- separati da una certa distanza $\mathbf{x}_{1,2}$ rispetto ad una distribuzione uniforme di materia in cui tutti gli oggetti sono "puntini" equidistanziati ("rumore poissoniano").
Tuttavia, il linguaggio nello spazio di Fourier è più diffuso, ed è il "frullatore" che amalgama i nostri ingredienti nella pratica quotidiana della ricerca in cosmologia. Il motivo fondamentale è che nello spazio di Fourier, i segnali complessi possono essere scomposti in somme di segnali semplici e fondamentali, con ciascun addendo indipendente dall'altro, così come il timbro di uno strumento musicale (il motivo per cui distinguiamo un flauto da un pianoforte anche quando suonano la stessa nota) può essere scomposto in sovrapposizione di frequenze. Questo porta ad una notevole semplificazione matematica dei calcoli, perché si può interpretare la dinamica delle perturbazioni dell'universo in questo modo:
Data una scala di lunghezza $L$ (solitamente misurata in Megaparsec, la dimensione caratteristica degli ammassi, circa trenta volte più grande delle dimensioni della via Lattea, il cui raggio è 30 kiloparsec), un modo nello spazio di Fourier è inversamente proporzionale alla scala ed è solitamente indicata con la lettera $k = \dfrac{2 \pi}{L} $. Modi piccoli corrispondono a grandi scale e strutture di grandi dimensioni, e viceversa. Le strutture si formano nei picchi di sovradensità, dove i potenziali gravitazionali sono più intensi, come illustrato dall'immagine qui sotto:
Lo stesso ragionamento si può fare per un'immagine in due dimensioni, pensate per esempio al volto di un essere umano:
Si può pensare che gli occhi, il naso e le orecchie siano "picchi" di una stessa onda di densità, così come le labbra sono invece l'effetto della sovrapposizione di più onde, che in parte si sommano facendo interferenza distruttiva, determinando solo una "struttura locale". Lo stesso ragionamento si può fare considerando i modi che si propagano in direzione verticale, dal mento alla fronte. Ecco che allora, raccogliendo tutti i modi che caratterizzano un'immagine $-$ per il cosmologo, la mappa tridimensionale della densità dell'universo $-$ si può lavorare nello spazio di Fourier per scomporre la complessità della ragnatela cosmica, e ricondurla poi alle sue proprietà statistiche. L'idea generale è "contare" l'abbondanza di modi per ciascun valore di scala, un po' come di un Paese del mondo possiamo raggruppare le città in un istogramma in base alla loro estensione od il numero di abitanti che hanno${^1}$, in modo più preciso, l'informazione minima che più ci interessa e la correlazione tra due punti qualunque della struttura a grande scala, cioè la statistica a due punti che chiamiamo spettro di potenza; non è l'unica, ma è la più diffusa e studiata, anche se c'è molta attenzione da parte dei ricercatori per estensioni e alternative. Più correlazione c'è, più ci si può aspettare di avere struttura organizzata.
$^{1}$ A rigore, quest'analogia è più appropriata per la funzione di massa degli aloni, non per lo spettro di potenza, che misura la varianza, cioè la deviazione dalla densità media, più che un valore puntuale di massa o estensione, scala per scala, ma credo che sia sufficiente per rendere l'idea, pur nella sua inaccuratezza.
I pacchetti $\texttt{halomod}$, $\texttt{hmf}$ consentono di generare lo "spettro di potenza" delle strutture cosmiche a partire dall'halo model. Questo modello assume che tutta la materia (materia oscura e barioni) sia contenuta dentro aloni di materia oscura (non esistono quindi né vuoti né fialmenti), e che le galassie a loro volta siano delle "palline" che occupano tali aloni in misura diversa a seconda che siano centrali nell'ammasso o satelliti, secondo modelli chiamati halo occupation distributions (HOD).
Rispetto allo spettro di potenza che si ottiene dalle equazioni della relatività generale (almeno al primo ordine perturbativo), l'halo model è automaticamente $\textbf{non lineare}$, ovvero tiene conto di come i processi di aggregazione delle galassie e degli ammassi a piccola scala (grandi k), determinino un eccesso di correlazione -- altrimenti detto, non è sufficiente descrivere una regione come sovrapposizione, cioè somma, di modi, occorre tenere conto anche del prodotto tra le varie onde di densità -- che si traduce in una "spalla" dello spettro di potenza a grande $k$.
Lo spettro dell'halo model si può definire in questo modo:
$$ P_{hm,totale}(k,z) = P_{hm, 1 halo}(k,z) + P_{hm, 2halo}(k,z)$$Dove:
$\phantom{a a a a}$
Si tenga conto che c'è poi un importante effetto non lineare, dovuto ai moti delle galassie in caduta nelle strutture cosmiche, chiamato redshift space distorsions (RSD) che l'halo model $\textbf{non}$ è in grado di riprodurre.
In generale, poi, ogni termine di 1halo e 2halo è un integrale sulla massa degli aloni ( = sommo i contributi degli aloni di massa, a partire dalla massa minima, che dipende dalle proprietà della materia oscura, più questa è "fredda", più sono i piccoli gli aloni che si possono formare). Questo perché la materia oscura "fredda" ha un'energia cinetica minore (di qui il suo nome) e tende a sfuggire meno ai potenziali gravitazionali, consentendo quindi la stabilizzazione di aloni potenzialmente fino a $10^{-6} M_{\odot}$ (un milionesimo della massa del sole, ovvero la massa della Terra!! Ma molto più rarefatta, un alone di quella massa avrebbe densità di $0.008 M_{\odot}/parsec^3$ e avrebbe un raggio di circa 2 mesi luce. Poiché il sistema solare, fino a Plutone, è grande 5 ore luce, questo vuol dire che l'alone sarebbe circa 240 volte più esteso del sistema solare). La somma si arresta attorno a $10^{18} M_{\odot}$, cioè gli aloni che formano gli ammassi di galassie; la massa tipica di una galassia, e dunque il suo alone, è $10^{11} - 10^{12} M_{\odot}$.
La formula di ciascun termine è la seguente:
$$ P_{\delta,\delta, 1h}(k,z) = \int_{M_{Min}} ^{M_{Max}} dM \ \frac{dN (M,z)}{dM} \left (\frac{\rho(k,M,z)}{\bar{\rho}(z)}\right)^2 \\ P_{\delta,\delta, 2h}(k,z) = P_{lin}(k,z) \left [\int_{M_{Min}} ^{M_{Max}} dM \ \frac{dN (M,z)}{dM} b(M,z) \left(\frac{\rho(k,M,z)}{\bar{\rho}(z)}\right) \right]^2$$Dove:
from IPython.display import Image
Image(filename='funzione_massa.png')
from IPython.display import Image
Image(filename='bias.png')
Ce ne sono di diversi tipi: il più semplice e fondamentale è il profilo Navarro-Frenk-White (NFW), che assume che gli aloni siano perfettamente sferici, con un tipico andamento di densità che cresce verso il centro dell'alone, ed una pendenza che cambia ad un raggio caratteristico funzione della massa $r_{s}$, che termina la sua "concentrazione". Più l'alone è di piccola massa, più piccolo è questo raggio, cioè è più "concentrato". Gli aloni a loro volta sono solitamente troncati ad un raggio massimo, il cosiddetto raggio viriale: il nome viene dal fatto che una struttura cosmica è stabile se è in equilibrio viriale, cioé la gravita ne determina l'aggregazione, ma non la fa collassare su se stessa.
from IPython.display import Image
Image(filename='profilo_densità.png')
from IPython.display import Image
Image(filename='profili_fourier.png')
from IPython.display import Image
Image(filename='profili_in_funzione_della_massa.png')
I due pacchetti $\texttt{halomod}$ e $\texttt{hmf}$ consentono di svogere tutti i calcoli necessari per estrarre gli spettri dall'halo model. Sotto inizializzeremo uno spettro con la ricetta "standard", che combina:
Genereremo una serie di spettri e simulazioni associate: il riferimento sarà il modello $\Lambda \rm{CDM}$, ma varieremo il parametro di densità per vedere gli effetti sulle simulazioni che genereremo successivamente.
# Cosmologia
import halomod
from halomod import DMHaloModel
import hmf
#HMF
from hmf import cosmology
from hmf import density_field
from hmf import Transfer
import powerbox as pbox
/home/andrea/anaconda3/envs/progetto_1/lib/python3.11/site-packages/halomod/halo_exclusion.py:18: UserWarning: Warning: Some Halo-Exclusion models have significant speedup when using Numba warnings.warn(
Spunti di riflessione: dipendenza dai parametri cosmologici, funzione di massa, profili degli aloni.
ref_Om_0 = hmf.cosmo.Planck15.Om0
# Generiamo cosmologia di riferimento
hm_ref = DMHaloModel()
redshift = 0.
hm_ref.z = redshift
hm_ref.cosmo_params['OmO'] = ref_Om_0
hm_ref.tracer_profile_model = "NFW" # Primo ingrediente: profilo di densità dell'alone di materia oscura
# Navarro - Frenk - White, lo standard più comune
hm_ref.bias_model = 'ST99' # Secondo ingrediente, il "bias dell'alone"
hm_ref.hmf_model = 'ST' # Terzo ingrediente: la funzione di massa, quanti (in numero) aloni abbiamo nell'universo per ogni valore di massa dell'alone?
modi_k = np.logspace(-4,3,200) # Generiamo più modi di quelli che halomod normalmente associa agli spettri (pur calcolando il valore dello spettro corrispondente), servirà dopo
# Genero lo spettro di potenza, termini 1halo, 2halo e complessivo
index_1h_max = np.where(hm_ref.power_1h_auto_matter_fnc(modi_k) == np.max(hm_ref.power_1h_auto_matter_fnc(modi_k)))[0]
max_val_1h = np.max(hm_ref.power_1h_auto_matter_fnc(modi_k)) # valore corrispondente
# Ridefinisco la funzione del 1halo a tratti, in modo che a grandi scale (piccoi k) restituisca un plateau
def pk_1halo_corretto(k,pk):
return np.piecewise(k, [k < modi_k[index_1h_max],k >= modi_k[index_1h_max] ], [lambda x: max_val_1h ,lambda x :pk.power_1h_auto_matter_fnc(x)])
# Generazione della simulazione di riferimento, LCDM
boxlength = 1000
N = 256
#print('Massimo k campionato: ',np.pi * N / boxlength, 'h / Mpc')
#print('Minima scala campionata:', boxlength/N, 'h / Mpc')
# Ovvero, scale più piccole delle scale tipiche degli ammassi di galassie, "vediamo" la non linearità
# dello spettro di potenza
# Sezione bidimensionale di una simulazione a dato redshift
lnpb_ref_tot = pbox.LogNormalPowerBox(
N=N, # Number of grid-points in the box
dim=2, # 2D box
pk = lambda k: hm_ref.power_auto_matter_fnc(k), # The power-spectrum
boxlength = boxlength , # Size of the box (sets the units of k in pk)
seed = 1010,ensure_physical=True
# Use the same seed as our powerbox
)
lnpb_ref_2h = pbox.LogNormalPowerBox(
N=N, # Number of grid-points in the box
dim=2, # 2D box
pk = lambda k: hm_ref.power_2h_auto_matter_fnc(k), # The power-spectrum
boxlength = boxlength , # Size of the box (sets the units of k in pk)
seed = 1010,ensure_physical=True
# Use the same seed as our powerbox
)
lnpb_ref_1h = pbox.LogNormalPowerBox(
N=N, # Number of grid-points in the box
dim=2, # 2D box
pk = lambda k: pk_1halo_corretto(k,hm_ref), # The power-spectrum
boxlength = boxlength , # Size of the box (sets the units of k in pk)
seed = 1010,ensure_physical=True
# Use the same seed as our powerbox
)
# Inizializziamo l'halo model
modi_k = np.logspace(-4,3,200) # Generiamo più modi di quelli che halomod normalmente associa agli spettri (pur calcolando il valore dello spettro corrispondente), servirà dopo
# Genero lo spettro di potenza, termini 1halo, 2halo e complessivo
index_1h_max = np.where(hm_ref.power_1h_auto_matter_fnc(modi_k) == np.max(hm_ref.power_1h_auto_matter_fnc(modi_k)))[0]
# indice del vettore con il massimo valore del termine di 1halo, quando raggiunge il plateau
# IL TERMINE 1 HALO NON FUNZIONA ESATTAMENTE, È UN PROBLEMA NUMERICO (SIGH, BENVENUTI NELLA QUOTIDIANITÀ DEL FISICO)!!!
plt.figure(figsize=(18,12))
plt.loglog(modi_k, hm_ref.power_auto_matter_fnc(modi_k), label='Spettro di potenza',lw=3)
plt.plot(modi_k, hm_ref.power_1h_auto_matter_fnc(modi_k), ls='--', label=' termine 1-halo ',lw=3)
plt.plot(modi_k, hm_ref.power_2h_auto_matter_fnc(modi_k), ls='--', label=' termine 2-halo ',lw=3)
plt.hlines(np.max(hm_ref.power_1h_auto_matter_fnc(modi_k)),modi_k[0],modi_k[index_1h_max] ,color = 'red', ls='--', label=' termine 1-halo come dovrebbe essere...',lw=3)
plt.grid()
plt.legend(fontsize=25)
plt.xlabel("k [$Mpc^{-1} h$]",size=32)
plt.ylabel(r" P(k)$ \ [{ Mpc^3}h^{-3}]$",size=32);
plt.ylim(1e-3,1e5)
/home/andrea/anaconda3/envs/progetto_1/lib/python3.11/site-packages/halomod/tools.py:653: UserWarning: You are using an un-normalized mass function and bias function pair.Bias ST99 has the following paired HMF model: (<class 'hmf.mass_function.fitting_functions.SMT'>,). Matter correlations are not well-defined. warnings.warn(
(0.001, 100000.0)
# Interludio: correggere le bizzarrie numeriche, al volo
max_val_1h = np.max(hm_ref.power_1h_auto_matter_fnc(modi_k)) # valore corrispondente
plt.figure(figsize=(18,12))
plt.plot(modi_k, pk_1halo_corretto(modi_k,hm_ref), ls='-', label=' termine 1-halo corretto ',lw=3)
plt.plot(modi_k, hm_ref.power_1h_auto_matter_fnc(modi_k), ls='--', label=' termine 1-halo ',lw=3)
plt.grid()
plt.legend(fontsize=25)
plt.xlabel("k [$Mpc^{-1} h$]",size=32)
plt.ylabel(r" P(k)$ \ [{ Mpc^3}h^{-3}]$",size=32);
plt.loglog()
plt.ylim(1e-4,1e4)
(0.0001, 10000.0)
Le simulazioni cosmologiche complete richiedono grande potenza di calcolo e la soluzione di complessi e lunghi calcoli per determinare la traiettoria delle particelle, sia che si modellino i soli effetti gravitazionali di particelle non collisionali (simulazioni N-Body, di fatto sola materia oscura), sia che si cerchi di tenere conto di effetti più complessi dati da particelle collisionali, capaci cioè di esercitare pressione, come i barioni, con le simulazioni idrodinamiche, che, in presenza di campi magnetici, onnipresenti nell'universo, diventano le ancor più complicate simulazioni magnetoidrodinamiche.
Tuttavia, si può accelerare il processo ricorrendo alle proprietà statistiche della struttura a grande scala! L'universo ha una densità media $\bar{\rho}$ che evolve nel tempo, rispetto alla quale le strutture e i vuoti rappresentano sovra e sotto densità, quantificate dal valore del contrasto di densità, ovvero la differenza tra la densità locale e quella media dell'universo che ci ospita.
Notare che la densità locale è funzione della posizione (per questo è locale) e del redshift, mentre la densità media è la stessa in tutto l'universo e cambia solo con il redshift $\delta(\textbf{x},z) = \dfrac{\rho(\textbf{x},z) - \bar{\rho}(z) }{\bar{\rho}(z)}$. Si tratta di un numero compreso tra - 1 (densità locale nulla) e, volendo, infinito, o arbitrariamente grande rispetto alla densità media (strutture molto molto dense concentrate). Ora, la cosmologia diventa non lineare quando tale contrasto diventa confrontabile con 1 o superiore, cosa che accade tipicamente:
Per avere una scala di riferimento, tenete presente che, ad oggi, il contrasto di densità vale
Notate che decresce man mano che si va a strutture più grandi, cioè, "zoomando" indietro, ci si avvicina sempre di più all'omogeneità con densità pari alla densità media del background (contrasto vicino a 0) ma, nel contempo, strutture ancora molto grandi sono di fatto lontane dall'ipotesi lineare.
Tale passaggio dalla linearità alla non linearità determina un cambiamento delle proprietà statistiche della struttura a grande scala. Originariamente, il valore dell'ampiezza di tali "semi" di contrasto di densità può essere considerato un numero casuale (notare non c'è un meccanismo deterministico conosciuto!) dato dalle fluttuazioni quantistiche dello spazio tempo di poco successive al Big Bang, è distribuito come una gaussiana ("normalmente"). La distribuzione normale del contrasto di densità $$f(x) = \frac{1}{\sqrt{ 2 \pi \sigma(x)^2}} e^{\frac{-x^2}{2\ \sigma(x)^2}} $$ ha media 0 (perché in media la densità è la densità media, quindi il contrasto medio è zero) e varianza $\sigma^2(x)$ funzione della posizione (o, nello spazio di Fourier, della scala $k$). Lo spettro di potenza ci dà in sostanza i valori della varianza ad ogni scala k e, siccome sappiamo che il contrasto in media è nullo, è la sola funzione che contiene tutte le informazioni necessarie per descrivere l'universo! Una singola funzione!
Ora, quando le strutture si formano, la distribuzione delle fluttuazioni cambia: in particolare, gli effetti della gravità a piccole scale $\textit{correlano}$.Questo vuol dire che la varianza non dipende solo dalla singola scala $k$, ma anche da eventuali coppie di scale $\sigma^2(k_1,k_2)$ (in un caso discreto, la matrice di covarianza diventa non più diagonale, lo spettro di potenza, che è la sua controparte continua, cambia forma e viene descritto dall'halo model). Si noti, en passant, che a grandi scale l'inflazione dopo il Big Bang ha un effetto simile, inducendo un livello di non-gaussianità primordiale non nullo.
Tale modificazione nelle proprietà statistiche dell'universo fa sì che la distribuzione della variabile casuale "contrasto di densità" si avvicini ad una $\textit{log}$-normale, $$f(x)=\frac{e^{-\frac{(\ln x-\mu)^2}{2\sigma^2}}}{x \sqrt{2\pi}{\sigma}}$$ cioè una distribuzione che rimane sempre compresa tra - 1 e infinito (là dove, tra l'altro, la distribuzione normale può assumere valori negativi arbitrariamente alti e non fisici, perché la varianza, cioè lo spettro di potenza, cresce al decrescere del redshift). Inoltre, siccome le strutture si formano prima a piccola scala, diventa automaticamente più probabile trovare eccessi di contrasto di densità, valori grandi e positivi rispetto ad una distribuzione normale. Mentre le code della distribuzione normale sono simmetriche, la lognormale ha lunghe code positive, catturando questo effetto, e nel contempo il sposta il proprio picco verso valori negativ del contrasto di densità. Si noti che la distribuzione lognormale resta un $\textbf{modello statistico}$, la vera distribuzione di materia potrebbe essere diversa e più complicata, al limite non disporre neppure di una forma analitica chiusa, di una "formula", ma già Hubble aveva notato che questa distribuzione descrive molto bene il posizionamento spaziale delle galassie, ed è ragionevolmente un buon modello empirico. Nell'immagine sotto vediamo bene come una distribuzione normale può venir progressivamente distorta in una gaussiana, e come le code crescano di lunghezza al crescere della varianza, mentre per varianza piccola la distribuzione è indistinguibile da una gaussiana simmetrica (la funzione è centrata in 1 perché qui viene mostrato il rapporto tra il valore della densità e la densità media :
Ora, si può ribaltare il discorso: se noi abbiamo un modello statistico (non fisico, mi raccomando) per la distribuzione della materia nell'universo, possiamo evitare di fare tutti i calcoli retrostanti una vera simulazione cosmologica e popolare un cubo di dati, che rappresenta una regione dell'universo, con dei valori di densità estratti da una distribuzione lognormale per la quale abbiamo specificato, ad un dato redshift, uno spettro di potenza. E, siccome si tratta di estrarre dei valori da una densità di probabilità, non risolvere equazioni, il processo è infinitamente più rapido, nell'ordine di qualche minuto: si può facilmente eseguire una simulazione su di un computer portatile.
D'altra parte, come sempre succede in fisica, quando le scale diventano molto non lineari, il meccanismo non funziona più: il contrasto di densità ha una frequenza di valori molto elevati e positivi maggiore di quanto predetto da una lognormale. Indicativamente, riusciamo a riprodurre la non-linearità in modo indistinguibile dalle simulazioni e dalla teoria fino a scale k ∼ $3 Mpc^{-1}$, cioè fino a circa la scala degli ammassi di galassie, con una lusinghiera precisione dell'1%.
boxlength = 1000
N = 256
print('Massimo k campionato: ',np.pi * N / boxlength, 'h / Mpc')
print('Minima scala campionata:', boxlength/N, 'h / Mpc')
# Ovvero, scale più piccole delle scale tipiche degli ammassi di galassie, "vediamo" la non linearità
# dello spettro di potenza
Massimo k campionato: 0.8042477193189871 h / Mpc Minima scala campionata: 3.90625 h / Mpc
In generale, una simulazione appare come una "scatola", una regione di cielo cubica o a forma di parallelepipedo.
Gli ingredienti necessari sono:
Inoltre, si può avere interesse a calibrare le dimensioni della simulazione in modo tale da campionare adeguatamente una data regione dello spettro, ad esempio le oscillazioni acustiche barioniche (BAO, con una dimensione tipica di 150 Mpc), ovvero riuscire a raccogliere sufficienti valori del campo di densità a quella scala affinché sia statisticamente robusto e affidabile (in altre parole: possiamo dire che tutte le mele sono rosse solo perché ne abbiamo vista una rossa? No, abbiamo bisogno di un numero ragionevole di campioni per trarre conclusioni sicure).
Inoltre, vogliamo "vedere nel dettaglio" le scale di interesse, risolverle con precisione in un'immagine rispetto alle altre scale, esattamente come la qualità di una macchina fotografica è determinata (anche) dalla sua risoluzione.
In generale, il più piccolo modo che può venir misurato è pari a $k_{\rm{min}} = 2\pi/L$ dove $L$ è la lunghezza del lato della scatola, mentre il modo più grande, pari alla frequenza di Nyquist del sistema, è pari a $k_{\rm{max}} = N_{\rm{griglia}}\pi/L$. Il modo $k_\rm{min}$ dà anche un'indicazione sulla risoluzione (dimensione del bin di campionamento), che devono essere un paio di volte almeno minore della scala che si vuole campionare con precisione (principio di Nyquist).
In altre parole:
Come spesso accade in fisica, ciò che si vuole fare è un $\textbf{compromesso}$ di tante variabili, e non c'è una risposta univoca, bensì $\textbf{dipende}$ da ciò che uno vuole fare!
Temi di discussione: lognormali, linearità, non linearità. Come si campiona?
# Sezione bidimensionale di una simulazione a dato redshift
plt.figure(figsize=(15,7))
plt.title('Complessivo',size=22)
plt.imshow(lnpb_ref_tot.delta_x(),extent=(0,boxlength,0,boxlength),cmap='viridis',alpha=1.,vmax=100)
plt.colorbar()
plt.show()
# Sezione bidimensionale di una simulazione a dato redshift
plt.figure(figsize=(15,7))
plt.title('2 halo',size=22)
plt.imshow(lnpb_ref_2h.delta_x(),extent=(0,boxlength,0,boxlength),cmap='viridis',alpha=1.,vmax=100)
plt.colorbar()
plt.show()
# Sezione bidimensionale di una simulazione a dato redshift
plt.figure(figsize=(15,7))
plt.title('1 halo',size=22)
plt.imshow(lnpb_ref_1h.delta_x(),extent=(0,boxlength,0,boxlength),cmap='viridis',alpha=1.,vmax=100)
plt.colorbar()
plt.show()
pb = pbox.LogNormalPowerBox(N=N, dim=3, pk = lambda k: hm_ref.power_auto_matter_fnc(k), boxlength= boxlength, seed=1234,ensure_physical=True)
deltax = pb.delta_x() # altra cosa importante: ci vuole del tempo...
plt.figure(figsize=(15,7))
plt.imshow(np.mean(deltax[:,:,:],axis=0),extent=(0,boxlength,0,boxlength),cmap='coolwarm',) # capire che cosa si intende
plt.colorbar()
plt.show()
# Estraggo lo spettro di potenza dalla simulazione,
p_k, kbins = pbox.get_power(deltax,pb.boxlength,)#,bins = 80)
plt.figure(figsize=(18,12))
plt.loglog(hm_ref.k_hm, hm_ref.power_auto_matter_fnc(hm_ref.k_hm), label='Analitico')
plt.plot(kbins,p_k,label="Campionato da simulazione")
plt.grid()
plt.legend(fontsize=25)
plt.xlabel("k [$Mpc^{-1} h$]",size=32)
plt.ylabel(r" P(k)$ \ [{ Mpc^3}h^{-3}]$",size=32);
plt.loglog()
#plt.ylim(1e-4,1e4)
#plt.xlim(2*1e-2,1)
[]
#hm_alt.cosmo
ref_Om_0 = hmf.cosmo.Planck15.Om0
ref_Ob_0 = hmf.cosmo.Planck15.Ob0
ref_H_0 = hmf.cosmo.Planck15.H0
# Generiamo cosmologia di riferimento
hm_alt = DMHaloModel(cosmo_params = {'Ob0':0.2,})
# 'Om0': 0.3,
# 'H0': 46.0,})
redshift = 0.
hm_alt.z = redshift
#hm_alt.cosmo_params['OmO'] = ref_Om_0#/2
#hm_alt.cosmo_params['ObO'] = ref_Ob_0#*15
#hm_alt.cosmo_params['H0'] = ref_H_0#*15
hm_alt.sigma_8 = 1.4
hm_alt.tracer_profile_model = "NFW" # Primo ingrediente: profilo di densità dell'alone di materia oscura
# Navarro - Frenk - White, lo standard più comune
hm_alt.bias_model = 'ST99' # Secondo ingrediente, il "bias dell'alone"
hm_alt.hmf_model = 'ST' # Terzo ingrediente: la funzione di massa, quanti (in numero) aloni abbiamo nell'universo per ogni valore di massa dell'alone?
modi_k = np.logspace(-4,3,200) # Generiamo più modi di quelli che halomod normalmente associa agli spettri (pur calcolando il valore dello spettro corrispondente), servirà dopo
# Genero lo spettro di potenza, termini 1halo, 2halo e complessivo
index_1h_max = np.where(hm_alt.power_1h_auto_matter_fnc(modi_k) == np.max(hm_alt.power_1h_auto_matter_fnc(modi_k)))[0]
max_val_1h = np.max(hm_alt.power_1h_auto_matter_fnc(modi_k)) # valore corrispondente
# Ridefinisco la funzione del 1halo a tratti, in modo che a grandi scale (piccoi k) restituisca un plateau
def pk_1halo_corretto(k,pk):
return np.piecewise(k, [k < modi_k[index_1h_max],k >= modi_k[index_1h_max] ], [lambda x: max_val_1h ,lambda x :pk.power_1h_auto_matter_fnc(x)])
# Generazione della simulazione di riferimento, LCDM
boxlength = 1000
N = 256
#print('Massimo k campionato: ',np.pi * N / boxlength, 'h / Mpc')
#print('Minima scala campionata:', boxlength/N, 'h / Mpc')
# Ovvero, scale più piccole delle scale tipiche degli ammassi di galassie, "vediamo" la non linearità
# dello spettro di potenza
# Sezione bidimensionale di una simulazione a dato redshift
lnpb_alt_tot = pbox.LogNormalPowerBox(
N=N, # Number of grid-points in the box
dim=2, # 2D box
pk = lambda k: hm_alt.power_auto_matter_fnc(k), # The power-spectrum
boxlength = boxlength , # Size of the box (sets the units of k in pk)
seed = 1010,ensure_physical=True
# Use the same seed as our powerbox
)
lnpb_alt_2h = pbox.LogNormalPowerBox(
N=N, # Number of grid-points in the box
dim=2, # 2D box
pk = lambda k: hm_alt.power_2h_auto_matter_fnc(k), # The power-spectrum
boxlength = boxlength , # Size of the box (sets the units of k in pk)
seed = 1010,ensure_physical=True
# Use the same seed as our powerbox
)
lnpb_alt_1h = pbox.LogNormalPowerBox(
N=N, # Number of grid-points in the box
dim=2, # 2D box
pk = lambda k: pk_1halo_corretto(k,hm_alt), # The power-spectrum
boxlength = boxlength , # Size of the box (sets the units of k in pk)
seed = 1010,ensure_physical=True
# Use the same seed as our powerbox
)
# Sezione bidimensionale di una simulazione a dato redshift
plt.figure(figsize=(15,7))
plt.title('Complessivo',size=22)
plt.imshow(lnpb_alt_tot.delta_x(),extent=(0,boxlength,0,boxlength),cmap='viridis',alpha=1.,vmax=100)
plt.colorbar()
plt.show()
/home/andrea/anaconda3/envs/progetto_1/lib/python3.11/site-packages/halomod/tools.py:653: UserWarning: You are using an un-normalized mass function and bias function pair.Bias ST99 has the following paired HMF model: (<class 'hmf.mass_function.fitting_functions.SMT'>,). Matter correlations are not well-defined. warnings.warn( /home/andrea/anaconda3/envs/progetto_1/lib/python3.11/site-packages/powerbox/powerbox.py:386: RuntimeWarning: invalid value encountered in log return np.log(1 + self.correlation_array())
# Sezione bidimensionale di una simulazione a dato redshift
plt.figure(figsize=(15,7))
plt.title('2 halo',size=22)
plt.imshow(lnpb_alt_2h.delta_x(),extent=(0,boxlength,0,boxlength),cmap='viridis',alpha=1.,vmax=100)
plt.colorbar()
plt.show()
# Sezione bidimensionale di una simulazione a dato redshift
plt.figure(figsize=(15,7))
plt.title('1 halo',size=22)
plt.imshow(lnpb_alt_1h.delta_x(),extent=(0,boxlength,0,boxlength),cmap='viridis',alpha=1.,vmax=100)
plt.colorbar()
plt.show()
pb_confronto = pbox.LogNormalPowerBox(N=N, dim=3, pk = lambda k: hm_alt.power_auto_matter_fnc(k), boxlength= boxlength, seed=1234,ensure_physical=True)
deltax_alt = pb_confronto.delta_x() # altra cosa importante: ci vuole del tempo...
plt.figure(figsize=(15,7))
plt.imshow(np.mean(deltax_alt[:,:,:],axis=0),extent=(0,boxlength,0,boxlength),cmap='coolwarm') # capire che cosa si intende
plt.colorbar()
plt.show()
# Estraggo lo spettro di potenza dalla simulazione,
p_k_alt, kbins = pbox.get_power(deltax_alt,pb.boxlength,)#,bins = 80)
plt.figure(figsize=(18,12))
plt.loglog(hm_ref.k_hm, hm_ref.power_auto_matter_fnc(hm_ref.k_hm), label='$\Lambda$CDM analitico')
plt.loglog(hm_alt.k_hm, hm_alt.power_auto_matter_fnc(hm_ref.k_hm), label='Cosmologia alternativa analitico')
plt.plot(kbins,p_k,label="$\Lambda$CDM, campionato da simulazione")
plt.plot(kbins,p_k_alt,label="Cosmologia alternativa, campionato da simulazione")
plt.grid()
plt.legend(fontsize=25)
plt.xlabel("k [$Mpc^{-1} h$]",size=32)
plt.ylabel(r" P(k)$ \ [{ Mpc^3}h^{-3}]$",size=32);
plt.loglog()
#plt.ylim(1e-4,1e4)
#plt.xlim(2*1e-2,1)
[]
# Vediamo se effettivamente abbiamo una distribuzione lognormale
plt.figure(figsize=(12,5))
bins = np.logspace(-7, np.log10(30.), 100)
_ = plt.hist((pb.delta_x()).flatten(), bins=bins,label='$\Lambda$CDM riferimento')
_ = plt.hist((pb_confronto.delta_x()).flatten(), bins=bins,label='cosmologia alternativa',alpha=0.8)
# format the axes
#plt.xscale('log')
plt.xlabel(r"$\delta$")
plt.ylabel(r"$N_\mathrm{cells}$")
plt.xlim(1e-4,100 )
plt.semilogx()
plt.legend()
<matplotlib.legend.Legend at 0x7fced05b1310>
Costruiamo ora una simulazione 3D per estrarre lo spettro di potenza
Temi di discussione: profondità in redshift, che senso ha uno spettro 3D a z fissato?
pb = pbox.LogNormalPowerBox(N=N, dim=3, pk = lambda k: hm_ref.power_auto_matter_fnc(k), boxlength= boxlength, seed=1234)
deltax = pb.delta_x() # altra cosa importante: ci vuole del tempo...
plt.figure(figsize=(15,7))
plt.imshow(np.mean(deltax[:,:,:],axis=0),extent=(0,boxlength,0,boxlength),cmap='coolwarm') # capire che cosa si intende
plt.colorbar()
plt.show()
# Estraggo lo spettro di potenza dalla simulazione,
p_k, kbins = pbox.get_power(deltax,pb.boxlength,)#,bins = 80)
plt.figure(figsize=(18,12))
plt.loglog(hm_ref.k_hm, hm_ref.power_auto_matter_fnc(hm_ref.k_hm), label='Analitico')
plt.plot(kbins,p_k,label="Campionato da simulazione")
plt.grid()
plt.legend(fontsize=25)
plt.xlabel("k [$Mpc^{-1} h$]",size=32)
plt.ylabel(r" P(k)$ \ [{ Mpc^3}h^{-3}]$",size=32);
plt.loglog()
#plt.ylim(1e-4,1e4)
#plt.xlim(2*1e-2,1)
[]
Sotto vediamo due grafici che rappresentano lo spettro di potenza delle strutture cosmiche: ha una bellezza quasi artistica, lo splendore della quantità di informazioni, la possibilità che ci dà di descrivere e riassumere la nostra (ancora incompleta, forse in parte erronea) comprensione dell'universo, mescolata all'eleganza della teoria retrostante, la relatività generale, racchiusa tuttavia in una grande semplicità. La versione a sinistra è molto significativa perché ci mostra come sia stato possibile effettuare varie misure, comprendo tantissime scale di interesse cosmologico ed astrofisico, ricorrendo a diverse campagne osservative e osservabili fisiche, dalla radiazione cosmica di fondo -- l'unica che finora ci ha consentito di misurare le scale più grandi, prima dell'inversione di pendenza -- alle misure a piccole scale ottenute con la foresta Lyman-$\alpha$.
A destra abbiamo una sua rappresentazione artistica, che ci mostra invece le sue proprietà salienti (la sua "topografia", caratterizzata dal punto di inversione e dal segnale ondulatorio delle oscillazioni acustiche barioniche - si noti che le BAO non sono oscillazioni perché il segnale che vediamo nello spazio di Fourier è ondulato!, ne parleremo meglio più avanti) e lo scompone negli andamenti diversi funzione della scala, che indagheremo in questa sezione.
import matplotlib.image as mpimg
# read images
img_A = mpimg.imread('aggiungimi.jpg')
img_B = mpimg.imread('spettro_artistico.png')
# display images
fig, ax = plt.subplots(1,2,figsize=(25,12))
fig.tight_layout()
ax[0].axis('off');
ax[1].axis('off');
ax[0].imshow(img_A,)
ax[1].imshow(img_B)
<matplotlib.image.AxesImage at 0x7fced02ac290>
Ma perché lo spettro di potenza ha la forma che vediamo? In generale, il suo andamento si può riassumere nella formula qui sotto:
$$ P(k,z) = D(z) \ T (k) \ P_{\rm{primordiale}} (k)$$Si tratta di una formulazione molto semplice e con la sua eleganza, in quanto ci sta dicendo che, fissato uno spettro primordiale, funzione solo del modo $k$, è possibile spostarci avanti e indietro nel tempo tramite la funzione di crescita $D(z)$, mentre la forma dello spettro ad ogni modo k viene modificata dalla funzione di crescita T(k). Si noti che le modifiche spaziali (fondamentalmente, cambio di pendenza dall'andamento rettilineo originario che ancora caratterizza le scale più grandi, ovvero i modi più piccoli) sono indipendenti dalle modifiche temporali, che fissano l'ampiezza dello spettro ad ogni redshift. Questa proprietà è un'altra manifestazione del concetto di linearità ma, sia ben chiaro, è vera appunto solo in regime lineare, dove è possibile disaccoppiare le interazioni tra le varie scale.
Ma perché lo spettro ha la forma che vediamo? Per capirlo, ricorreremo ad un altro strumento fondamentale del cosmologo, il Boltzmann solver $\texttt{CAMB}$. Un Boltzmann solver è un codice che combina l'evoluzione delle perturbazioni nello spazio (di Fourier) e nel tempo, dettata dalla relatività generale, con l'equazione di Boltzmann. Questa equazione ha come incognita una funzione di di distribuzione, solitamente denotata con $f(\textbf{x},z)$ o $n(\textbf{x},z)$, che può essere interpretata come una densità numerica (particelle per unità di volume). Per ogni specie di materia che popola l'universo, vogliamo la sua distribuzione di densità ad ogni dato istante e posizione, perché da questa è possibile ricavare la densità e la pressione del fluido, accoppiate tramite un'equazione di stato che consente poi di risolvere le equazioni per le perturbazioni di densità.
Sembra facile? Non lo è per nulla, perché la funzione di distribuzione dipende da energia e temperatura del fluido, da interazioni che inducono annichilazioni e decadimenti, dunque creano e distruggono particelle, specialmente nell'universo immediatamente successivo al Big Bang, a loro voltaf funzioni dell'espansione dell'universo (effetti di background) e dalle interazioni nelle strutture (effetti perturbativi). Il Boltzmann solver risolve quindi in modo accoppiato, ad ogni redshift, perturbazioni e equazione di Boltzmann, consentendoci di stabilire lo stato dinamico e termodinamico dell'universo ad ogni istante, con una precisione che dipende da alcune condizioni, fondamentalmente collegate, come sempre, alla distinzione tra regimi lineari e non lineari.
Per quanto complicata, l'evoluzione dei fluidi cosmologici si può però comprendere in modo relativamente semplice introducendo il concetto di disaccoppiamento (decoupling); ovvero, le varie specie (barioni, fotoni, materia oscura, neutrini...) interagiscono tra di loro fino a quando sono sufficientemente densi e dunque le particelle hanno un'elevata probabilità di scontrarsi, scambiandosi energia e quantità di moto. La probabilità di scontro dipende dal tipo di particella, fotoni e barioni sono molto interagenti, la materia oscura o i neutrini molto poco; pertanto, quanto più l'universo si espande, tanto più aumenterà il loro cammino libero medio, la distanza media che possono percorrere.
Ad un certo redshift, una data specie cessa di interagire (se non per effetto della gravità) e si disaccoppia; da quel momento, la sua temperatura scalerà soltanto con il redshift, $T(z) = \frac{T_\rm{disacc}}{1+z}$, perché non ci saranno più interazioni che possono modificarla, a parte l'espansione dell'universo che ne "diluisce" la temperatura. Al momento, tutte le specie cosmiche originarie sono di fatto disaccoppiate, tranne barioni e materia oscura, che, assieme, costituiscono la ragnatela cosmica grazie agli effetti della gravità.
Il fenomeno più importante nella storia primordiale del cosmo, a redshift circa 1100, ovvero circa 300 mila anni dopo il Big Bang, è il disaccoppiamento tra barioni e fotoni, che determina la formazione della radiazione cosmica di fondo (fotoni primordiali disaccoppiati) e fissa nello spettro di potenza le oscillazioni acustiche barioniche. Altri effetti di interazione tra materia oscura, barioni e "free streaming" dei fotoni, cioè il loro "sfuggire" dai potenziali gravitazionali, progressivamente "spianandoli", determina il cosiddetto "processamento" dello spettro di potenza descritto dalla funzione di trasferimento $T(k)$.
Esploreremo alcuni di questi effetti nelle prossime celle, partendo dal concetto di "entrata nell'orizzonte" delle varie scale spaziali.
Le varie scale spaziali hanno una dimensione che può essere più grande o più piccola della scala dell'orizzonte di Hubble -- qualitativamente, la scala più grande che, a ogni dato redshift, i fotoni possono raggiungere, ora attorno ai 4 Gigaparsec-- portando così informazione e dunque connessione causale, in particolare, interazione gravitazionale.
Le scale più grandi dell'orizzonte sono "congelate", cioè i potenziali gravitazionali non variano nel tempo, e le perturbazioni non possono collassare, non possono, cioè, determinare la formazione di strutture. Con il passare del tempo, sempre più scale diventano più piccole dell'orizzonte: entrano cioè in connessione causale, e possono cominciare ad interagire. Questo determina progressivamente la formazione delle strutture, prima le scale piccole poi le scale più grandi (si parla di teoria gerarchica, o clustering gerarchico, analogo all'evoluzione della funzione di massa al decrescere del redshift), purché queste siano più grandi della scala di Jeans, una scala minima al di sotto della quale (ovvero, per strutture più piccole di $L_{\rm{Jeans}}$) una struttura non può formarsi perché la pressione dei barioni si oppone con successo al collasso, disperdendo la perturbazione.
Un evento fondamentale nella storia dell'universo è il passaggio da un universo primordiale dominato dalla radiazione (quando, cioè, la densità dei fotoni è maggiore della densità della materia barionica e oscura) a quella dominata dalla materia (quando la densità dei fotoni diventa più piccola della densità di materia barionica ed oscura). Questo determina un fondamentale cambiamento nella crescita delle perturbazioni, che diventa più lento, determinando la formazione del punto di inversione dello spettro di potenza stesso. (Si noti che, da qualche miliardo di anni,siamo entrati nella fase in cui è l'energia oscura a dominare; l'orizzonte cioè progressivamente si ridurrà, determinando nuovamente un'uscita delle scale dall'orizzonte cosmologico.)
Vedremo l'evoluzione di questi effetti, scala per scala e fluido per fluido, nelle prossime celle.
# Prepariamo CAMB
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import sys, platform, os
import matplotlib
from matplotlib import pyplot as plt
import numpy as np
import camb
from camb import model, initialpower
# Latex: per impostare la grafica "come se fosse un vero articolo scientifico"
font = {'size' : 18, 'family':'STIXGeneral'}
axislabelfontsize='x-large'
plt.rc('font', **font)
plt.rcParams['text.usetex'] = True
np.set_printoptions(precision=3,linewidth=200,suppress=False)
# Inizializziamo la calcolatrice, come abbiamo già visto sopra
pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122) # fissiamo i parametri fondamentali della cosmologia: parametro
# di Hubble, densità dei barioni, densità della materia oscura
pars.InitPower.set_params(ns=0.965) # indice dello spettro primordiale
pars.set_matter_power(redshifts=[0., 0.8], kmax=2.0) # fissiamo la scala e i redshift ai quali calcolare gli spettri
# Spettri lineari
pars.NonLinear = model.NonLinear_none
results = camb.get_results(pars)
kh, z, pk = results.get_matter_power_spectrum(minkh=1e-4, maxkh=1, npoints = 200)
s8 = np.array(results.get_sigma8())
# Spettri non lineari: calcolati a partire da dei modelli estratti dalle simulazioni, chiamati "Halofit", cioè
# fit degli aloni di materia oscura individuati nelle simulazioni
pars.NonLinear = model.NonLinear_both
results.calc_power_spectra(pars)
kh_nonlin, z_nonlin, pk_nonlin = results.get_matter_power_spectrum(minkh=1e-4, maxkh=1, npoints = 200)
Note: redshifts have been re-sorted (earliest first)
data = camb.get_background(pars)
# Calcoliamo effettivamente le varie scale all'evolvere del redshift rispetto all'orizzonte di Hubble: a=1 è
# z= 0, cioè "oggi"
z = np.linspace(0,5000,5000)
ks = [1e-4,1e-2,1e-1, 1]
ev = data.get_redshift_evolution(ks, z, ['delta_baryon','delta_cdm', 'delta_photon','H','a','k/h'])
_, axs= plt.subplots(1,len(ks), figsize=(18,5))
for i, ax in enumerate(axs):
ax.loglog(ev[1,:, -2],1/ev[i,:, -3]*ev[1,:, -2],lw=2,)
ax.loglog(ev[1,:, -2],2*np.pi/(ev[i,:, -1]/pars.h)*ev[1,:, -2],ls='--',lw=2,)
ax.set_title(r'$k= %s \ \rm{Mpc^{-1}}$'%ks[i])
ax.set_xlabel('$ a $');
ax.axvline(1/(1+3300),0,1e4,ls='-.',color='red',lw=2,)
ax.grid()
axs[0].set_ylabel(r'scala \rm{[Mpc]}');
plt.legend(['scala fisica','scala dell\'\ orizzonte di Hubble','$a$ equivalenza materia-radiazione'], ncols=3,bbox_to_anchor=(-0.2,-0.15),);
/home/andrea/anaconda3/envs/progetto_1/lib/python3.11/site-packages/matplotlib/scale.py:253: RuntimeWarning: overflow encountered in power return np.power(self.base, a)
Al crescere di di $k$ (decrescere della scala di lunghezza fisica), le scale "entrano (ovvero, la retta arancione interseca quella blu) nell'orizzonte di Hubble" prima ( $a$ più piccolo, ovvero redshift $z$ più alto, perché $a = \frac{1}{1+z}$. Quando le scale entrano nell'orizzonte, sono in connessione causale, la gravità può "lavorare" e farle interagire. Notare il cambio di pendenza dell'orizzonte di Hubble. In particolare, distinguiamo tra le scale che attraversano l'orizzonte prima o dopo il redshift dell'equivalenza ($z \sim 3300$), denotato dalla retta verticale rossa
z = np.linspace(3000,4000,1000)
ks = [1e-4,1e-2,0.08,1e-1, 0.8,1]#0.001,0.1,1,100,0.01]
prefattore = 2.1 * 1e-9/np.array(ks)**3 * 2 * np.pi**2
ev = data.get_redshift_evolution(ks, z, ['delta_baryon','delta_cdm', 'delta_photon','delta_tot'])
_, axs= plt.subplots(1,len(ks), figsize=(18,5))
#print(ev)
for i, ax in enumerate(axs):
ax.loglog(z,ev[i,:, 0]**2*prefattore[i])
ax.loglog(z,ev[i,:, 1]**2*prefattore[i])
ax.loglog(z,ev[i,:, 2]**2*prefattore[i])
ax.loglog(z,ev[i,:, 3]**2*prefattore[i])
ax.axvline(1+3300,0,1e4,ls='-.',color='red',lw=2,)
ax.grid()
ax.set_title(r'$k= %s \ \rm{Mpc^{-1}}$'%ks[i])
ax.set_xlabel('$z$');
axs[0].set_ylabel('$\delta$');
plt.legend([r'$\Delta_{barioni}$', r'$\Delta_{{materia \ oscura}}$', r'$\Delta_\gamma$',r'$\Delta_{tot}$','$z$ equivalenza materia-radiazione'], ncols=5,bbox_to_anchor=(-0.2,-0.15),);
plt.semilogx()
[]
Ad alto redshift, prima dell'equivalenza radiazione - materia, gli effetti si distinguono in questo modo:
z = np.linspace(0,2000,1000)
ks = [1e-4,1e-2,0.08,1e-1, 0.8,1]#0.001,0.1,1,100,0.01]
ev = data.get_redshift_evolution(ks, z, ['delta_baryon','delta_cdm', 'delta_photon','delta_tot'])
_, axs= plt.subplots(1,len(ks), figsize=(25,8))
#print(ev)
for i, ax in enumerate(axs):
ax.loglog(z,ev[i,:, 0]**2*prefattore[i])
ax.loglog(z,ev[i,:, 1]**2*prefattore[i])
ax.loglog(z,ev[i,:, 2]**2*prefattore[i])
ax.loglog(z,ev[i,:, 3]**2*prefattore[i])
ax.axvline(1+3300,0,1e4,ls='-.',color='red',lw=2,)
ax.axvline(1+1100,0,1e4,ls='-.',color='midnightblue',lw=2,)
ax.grid()
ax.set_title(r'$k= %s \ \rm{Mpc^{-1}}$'%ks[i])
ax.set_xlabel('$z$');
axs[0].set_ylabel('$\delta$');
plt.legend([r'$\Delta_{barioni}$', r'$\Delta_{{materia \ oscura}}$', r'$\Delta_\gamma$',r'$\Delta_{tot}$','$z$ equivalenza materia-radiazione','$z$ disaccoppiamento radiazione - materia'], ncols=6,bbox_to_anchor=(-0.2,-0.15),);
plt.semilogx()
[]
A basso redshift, dopo l'equivalenza radiazione - materia, gli effetti si distinguono in questo modo:
L'andamento delle perturbazioni in redshift si traduce nell'andamento per scala riassunto dallo spettro di potenza. Man mano che il tempo passa, più scale entrano nell'orizzonte e lo spettro di potenza viene smorzato da destra (piccole scale) a sinistra, con la contemporanea progressiva attivazione delle oscillazioni acustiche.
redshift = np.array([0,1e-4,100,500, 1000,3300,1e4,1e5])#np.array([4000, 3000, 1000, 500, 0])
scale = np.logspace(-4,2,100)
prefattore = 2.1 * 1e-9/scale**3 * 2 * np.pi**2
pk_primordiale = results.Params.scalar_power(scale)
indice_p_max = np.zeros(len(redshift))
plt.figure(figsize=(12,5))
plt.title('Materia')
for j, redshifts in enumerate(redshift):
spettro = np.array(data.get_redshift_evolution(scale, redshifts, ['delta_tot']).flatten()**2)*prefattore
plt.loglog(scale,spettro, label = r'$z= %s $'%int(redshifts) )
#print(np.where(spettro == np.max(spettro)))
indice_p_max[j] = int(np.where(spettro == np.max(spettro))[0])
plt.ylabel(r'$k \ [Mpc^{-1} \ h ]$')
plt.ylabel(r'$\delta^2$')
plt.loglog(scale,pk_primordiale * scale ** 4 / (scale**3/(2*np.pi**2)),label='primordiale')
plt.legend(ncols=1,bbox_to_anchor=(1,1),)
plt.grid()
plt.title('Scala del picco ')
plt.scatter(redshift,[scale[int(index)] for index in indice_p_max])
plt.axhline(scale[int(indice_p_max[0])],0,np.max(redshift),ls='-',color='red')
plt.axvline(1+3300,0,1e4,ls='-.',color='black',lw=2,)
plt.semilogx()
plt.grid()
plt.xlabel(r'$z$')
plt.ylabel(r'$k \ [Mpc^{-1} \ h ]$')
Text(0, 0.5, '$k \\ [Mpc^{-1} \\ h ]$')
Il picco dello spettro si sposta verso scale più grandi (verso sinistra) fino all'equivalenza radiazione materia o per redshift poco più bassi (nessun effetto fisico è davvero istantaneo!), poi si stabilizza perché le perturbazioni crescono tutte assieme, senza differenza tra l'essere sopra o sotto l'orizzonte.
redshift = np.array([0,1e-4,100,500, 1000,3300,1e4,1e5])#np.array([4000, 3000, 1000, 500, 0])
scale = np.logspace(-4,2,100)
prefattore = 2.1 * 1e-9/scale**3 * 2 * np.pi**2
pk_primordiale = results.Params.scalar_power(scale)
indice_p_max = np.zeros(len(redshift))
plt.figure(figsize=(12,5))
plt.title('Radiazione (fotoni)')
for j, redshifts in enumerate(redshift):
spettro = np.array(data.get_redshift_evolution(scale, redshifts, ['delta_photon']).flatten()**2)*prefattore
plt.loglog(scale,spettro, label = r'$z= %s $'%int(redshifts) )
indice_p_max[j] = int(np.where(spettro == np.max(spettro))[0])
plt.loglog(scale,pk_primordiale * scale ** 4 / (scale**3/(2*np.pi**2)),label='primordiale')
plt.legend(ncols=1,bbox_to_anchor=(1,1),)
plt.xlabel(r'$z$')
plt.ylabel(r'$\delta^2$')
plt.grid()
La fenomenologia di un fluido relativistico è ben diversa: le perturbazioni crescono al decrescere del redshift, ma i fotoni non formano struttura. Quando le distanze sono diventate troppo grandi perché fotoni e barioni possano interagire ( $z\sim 1100$, epoca del disaccoppiamento e della formazione della CMB) i fotoni sono liberi di muoversi per tutto l'universo, "sparpagliandosi" in modo praticamente uniforme (free streaming). Le oscillazioni vengono a loro volta progressivamente smorzate e ogni forma di correlazione (=struttura) viene persa (i valori in ordinata dello spettro di potenza sono ben più piccoli dello spettro di materia).
Le prossime due celle mostrano un video, anzi, un film, un film tanto semplice quanto può essere la modellizzazione fisica -- che riesce a ricondurre fenomeni complessi ad una quantità ridotta ('riduzionismo') di concetti elementari (potenziali, velocità, densità degli elementi del sistema...) -- quanto incredibilmente ricco e significativo, a mio parere di folgorante bellezza: la storia dell'universo, delle strutture che si formano, al passare del tempo, "narrata" a partire dall'evoluzione dei potenziali gravitazionali.
Può apparire incredibilmente semplice, ma le forze concorrenti nell'universo sono fondamentalmente di tre categorie soltanto:
# Per quello che interessa a noi, tutto quello che segue è solo codice preparatorio, senza un vero interesse
# per il significato delle varie operazioni svolte.
data= camb.get_results(pars)
zstar=data.get_derived_params()['zstar']
etastar=data.conformal_time(zstar)
eta = np.logspace(-1.7,np.log10(etastar),400)
k=10**np.linspace(-3,1,1200)
z = data.redshift_at_conformal_time(eta)
back_ev = data.get_background_time_evolution(eta, ['x_e', 'visibility'])
x_e=back_ev['x_e']
vis=back_ev['visibility']
adotoa = data.h_of_z(z)/(1+z) #
rho=data.get_background_densities(1/(1+z))
Rmat=(rho['baryon']+rho['cdm'])/(rho['photon']+rho['neutrino']+rho['nu'])
grhonu=data.grhormass[0]+data.grhornomass
Rv=grhonu/(grhonu+data.grhog)
omega = (data.grhob+data.grhoc)/np.sqrt(3*(data.grhog+grhonu))
Psi_init = -10/(4*Rv+15)
import camb.symbolic as cs
ev=data.get_time_evolution(k, eta, [cs.Psi_N,
'delta_photon',
cs.make_frame_invariant(cs.Delta_g, 'Newtonian'),
cs.make_frame_invariant(cs.v_b, 'Newtonian')
], lAccuracyBoost=4)
Psi = ev[:,:,0]
Delta_g = ev[:,:,1]
Delta_g_N=ev[:,:,2]
v_b_N = ev[:,:,3]
for i, etai in enumerate(eta):
for j, kj in enumerate(k):
if etai*kj < 0.1 and etai*omega < 0.1:
Psi[j,i]=Psi_init - 25/8*omega*etai *(8*Rv-3)/(4*Rv+15)/(2*Rv+15)
Delta_g_N[j,i] = Delta_g[j,i] - 4*adotoa[i]*( Psi_init/2*etai - 15*omega*etai**2/8*(4*Rv-5)/(4*Rv+15)/(2*Rv+15))
#pip install ffmpeg-python # decommentami su un altro computer per installare
# Anche qui, preparazione e plot, a noi interessa soprattutto il video,
import math
from matplotlib import animation
def latex_sci_format(num, dec=3):
return "{:.{dec}f}\\times 10^{{{}}}".format(num / (10 ** int(math.floor(math.log10(abs(num))))), int(math.floor(math.log10(abs(num)))), dec=dec)
fig, ax = plt.subplots()
plot_vars = {r'$\Psi$':Psi,
r'$\frac{1}{4}\Delta^{\rm Newt}_\gamma +\Psi$': Delta_g_N/4+Psi,
r'$\frac{1}{\sqrt{3}}v_b^{\rm Newt}$': v_b_N/np.sqrt(3) }
anim_fps=8
lines=[]
for lab in plot_vars.keys():
line, = ax.semilogx([], [], label=lab)
lines.append(line)
ax.axhline(0,color='k')
plt.legend(loc='upper left')
ax.set_xlim(k[0], k[-1])
ax.set_ylim(-1, 1)
ax.set_xlabel(r'$k\, {\rm Mpc}$')
# Define the update function
def update(i):
ax.set_title( r'$z= %s, \eta/{\rm Mpc}= %.2f, x_e = %.2f, \rho_m/\rho_{\rm rad}=%.2f$'
% (latex_sci_format(z[i]), eta[i], x_e[i], Rmat[i]))
for line, var in zip(lines,plot_vars.values()):
line.set_data(k, var[:, i])
return lines
# Create the animation
transfer_anim =animation.FuncAnimation(fig, update, frames=range(len(eta)), blit=True)
#writer = animation.writers['ffmpeg'](fps=anim_fps, bitrate=-1)
# can save like this
#transfer_anim.save(r'Newtonian_transfer_evolve_monopole_phi_vb.mp4', writer=writer, dpi=240)
See the output mp4 file here or play below:
Il video sopra è molto denso: in generale, mostra il comportamento dei potenziali gravitazionali ($\Psi$), delle perturbazioni dei fotoni ($\frac{1}{4}\Delta^{\rm Newt}_\gamma +\Psi$'), che interagiscono con i barioni, e delle velocità in funzione del redshift (che decresce, cioè parte da un universo molto primordiale e arriva fino alla formazione della CMB a z $\sim$ 1100) e delle velocità dei barioni $\frac{1}{\sqrt{3}}v_b^{\rm Newt}$.
Si noti:
# Si può anche eseguire a parte in questa cella
from IPython.display import display, Video
video = transfer_anim.to_html5_video()
display(Video(data=video, embed=True, width = 600))
# Il codice sotto si concentra sui soli potenziali gravitazionali e i valori del contrasto di densità, e genera nuovamente
# un video.
eta=np.logspace(-1.7,np.log10(data.tau0),500)[:-1]
k=10**np.linspace(-3,2,2000)
z = data.redshift_at_conformal_time(eta)
import camb.symbolic as cs
ev=data.get_time_evolution(k, eta, [cs.Psi_N,
'delta_cdm', # sychronous gauge cdm perturbation
'delta_photon',
'delta_baryon'], lAccuracyBoost=4)
Psi = ev[:,:,0]
Delta_c = ev[:,:,1]
Delta_g=ev[:,:,2]
Delta_b=ev[:,:,3]
x_e = data.get_background_time_evolution(eta, ['x_e'])['x_e']
adotoa = data.h_of_z(z)/(1+z) # comoving Hubble
rho=data.get_background_densities(1/(1+z))
Rmat=(rho['baryon']+rho['cdm'])/(rho['photon']+rho['neutrino']+rho['nu'])
zstar=data.get_derived_params()['zstar']
# series for super-horizon Psi
grhonu=data.grhormass[0]+data.grhornomass
Rv=grhonu/(grhonu+data.grhog)
omega = (data.grhob+data.grhoc)/np.sqrt(3*(data.grhog+grhonu))
# initial value of super-horizon Psi for unit curvature and no isocurvature
Psi_init = -10/(4*Rv+15)
for i, etai in enumerate(eta):
for j, kj in enumerate(k):
if etai*kj < 0.1 and etai*omega < 0.1:
Psi[j,i]=Psi_init - 25/8*omega*etai *(8*Rv-3)/(4*Rv+15)/(2*Rv+15)
fig, ax = plt.subplots()
plot_vars = {r'$\Delta_c$': Delta_c, r'$\Delta_b$': Delta_b,
r'$\Delta_\gamma$': Delta_g,r'$\Psi$': Psi}
anim_fps=8
lines=[]
for lab in plot_vars.keys():
line, = ax.loglog([], [], label=lab)
lines.append(line)
plt.axhline(0, ls='-', color='k')
plt.legend(loc='upper left')
ax.set_xlim(k[0], k[-1])
ax.set_xlabel(r'$k\, {\rm Mpc}$')
ax.set_ylim(1e-4, 2e5)
lines.append(ax.loglog([], [],ls='--', color='C1')[0])
lines.append(ax.loglog([], [],ls='--', color='C2')[0])
lines.append(ax.loglog([], [],ls='--', color='C3')[0])
# Define the update function
def update(i):
ax.set_title( r'$z= %s, \eta/{\rm Mpc}= %.2f, x_e = %.2f, \rho_m/\rho_{\rm rad}=%.2f$'
% (latex_sci_format(z[i],2), eta[i], x_e[i], Rmat[i]), fontsize=10)
for line, var in zip(lines,plot_vars.values()):
line.set_data(k, var[:, i])
lines[-3].set_data(k,-Delta_b[:, i])
lines[-2].set_data(k,-Delta_g[:, i])
lines[-1].set_data(k,-Psi[:, i])
if z[i] < zstar: # hide photons after recombination as irrelevant for matter
lines[2].set_data([],[])
lines[-2].set_data([],[])
return lines
# Create the animation
transfer_anim =animation.FuncAnimation(fig, update, frames=range(len(eta)), blit=True, repeat = False)
writer = animation.writers['ffmpeg'](fps=anim_fps, bitrate=-1)
transfer_anim.save(r'matter_transfer_evolve_sync.mp4', writer=writer, dpi=240)
/tmp/ipykernel_4984/1669222679.py:12: UserWarning: Data has no positive values, and therefore cannot be log-scaled. ax.set_xlim(k[0], k[-1])
This the video produced of the synchronous transfer function evolution (dashed lines are negative):
Il video può essere visto andando all'ultima cella di questo sito https://camb.readthedocs.io/en/latest/CAMBdemo.html
Il video è molto simile al precedente: più nel dettaglio, vediamo i valori di contrasto di densità $\Delta$ delle varie specie di materia e radiazione. Si noti che:
Temi: BAO, misura di $H_{0}$, righelli e candele standard, tensioni, inferenza statistica, cross-correlazioni, bias e traccianti.
Lo spettro di potenza non ci racconta soltanto la storia nel tempo e la struttura nello spazio dell'universo: è anche un fondamentale strumento di misura, che consente di ricavare una serie di quantità fondamentali per il modello $\Lambda$CDM: in particolare, la varianza delle fluttuazioni entro un raggio di 8 Mpc, $\sigma_8$, il parametro di densità della materia $\Omega_{m}$ dalla scala del picco (o punto di inversione, o turnover) e, dalle oscillazioni acustiche barioniche, misurare indirettamente il valore della costante di Hubble $H_{0}$. Il fatto che sia le BAO sia il punto di inversione siano legati a distanze invarianti è la ragione del nome che viene loro solitamente dato di "righelli standard". Sono indicatori di distanza a grandissima scala che si comportano in modo analogo alle "candele standard" dell'universo locale (Cefeidi, Supernove Ia...), oggetti la cui luminosità intrinseca è nota e di cui, data la luminosità apparente misurata a terra, è possibile calcolare la distanza e calibrare la scala delle distanze cosmologiche.
La scala tipica delle oscillazioni acustiche barioniche, 150 Mpc, è il più grande righello standard a nostra disposizione al momento del quale disponiamo da una quindicina d'anni di misure. La misura di $H_0$ che si ottiene non è però indipendente: è bensì funzione (si dice "degenere") del valore dell'orizzonte acustico $r_s$ che viene misurato con la radiazione cosmica di fondo, che va quindi fissato indipendentemente. In alternativa, si può misurare il rapporto tra la loro dimensione radiale ed angolare, ricorrendo alle candele standard nell'universo "locale" (basso redshift, $z$ compreso tra 0 e 2 da contrastare con $z=1100$ della CMB). Questo impedisce però di misurare direttamente $H_{0}$; per ritrovarla, dobbiamo imporre un modello di cosmologia, con i parametri nuovamente fissati dalla CMB.
Quello che risulta però è che le due misure di $H_0$ sono ad oggi in tensione: la CMB ci fornisce un valore (ottenuto direttamente, ma a partire da un'osservabile "complicata" come lo spettro della radiazione cosmica di fondo) più basso, $67.66\pm0.42$ $\rm{km \ s ^{-1} \ Mpc^{-1}}$, rispetto al valore di $73.04\pm1.04$ $\rm{km \ s ^{-1} \ Mpc^{-1}}$ ottenuto indirettamente, ma a partire da più semplici misure di distanza di oggetti astrofisici.
Una simile tensione, ma più debole, si osserva anche sulla normalizzazione dello spettro di potenza $\sigma_8$, che risulta invece più grande per la CMB e più piccola per lo spettro di potenza misurato con la struttura a grande scala.
Ma che cosa vuol dire che due misure fisiche sono "in tensione"? Per capirlo, tenteremo per prima cosa il fit dello spettro di potenza ed una misura indiretta delle oscillazioni acustiche barioniche, con il cosiddetto formalismo di Alcock-Paczinski.
L'idea è la seguente: se la nostra misura di distanza, dipendente dal modello cosmologico "fiduciale", è affidabile, allora la dipendenza dello spettro di potenza da $k$ può essere parametrizzata in questo modo: $P(k) = P(k/\alpha)$. I dati delle misure devono disporsi allora sulla curva in modo tale che tale curva riproduca lo spettro di potenza teorico; questo è vero solo se $\alpha = 1$. Ma questo vuol dire, indirettamente, che il valore di $H_{0}$ che stiamo assumendo non è corretto.
Di seguito, non ricaveremo $H_{0}$, che richiede tempo ed è complesso, ci limiteremo soltanto a trovare il valore di $\alpha$. Vedremo però che ogni dato fisico è affetto da incertezza, e questa incertezza rende la situazione più sfumata: non si può dire in assoluto che una misura sia "giusta" o "sbagliata", va confrontata con le incertezze, e questo implica il confronto tra modelli, ed il nostro grado di fiducia in essi, così come la nostra necessità di misure sempre più precise ed accurate, per scartare quante più alternative possibili.
Questo tipo di ragionamento si chiama inferenza statistica, ed è fondamentale in ogni branca della fisica: è il nostro unico criterio di verità, che combina matematica, statistica, precisione e buon senso, tutti strumenti fondamentali per fare il fisico!
# Visualizzazione fit dei dati
from getdist import plots, MCSamples
import getdist
k_data = kbins[kbins <= 1] #taglio i dati rumorosi a grandi scale, che comunque non servono
cut_index_in = 0 ##np.where((hm.k_hm >= 2*1e-2 ) & (hm.k_hm <= 1))[0][0]
cut_index_fin = np.where(kbins <= 1)[0][-1]
# Il nostro modello di spettro ha due componenti: un termine con le BAO e uno senza.
spettro_BAO = hm_ref # quello con le BAO è l'halo model...
hm_ref.z
0.0
# Quello senza BAO è costruito con l'halo model, ma cambiamo la funzione di trasferimento.
spettro_noBAO = DMHaloModel()
spettro_noBAO.transfer_model = 'EH_NoBAO'
spettro_noBAO.tracer_profile_model = "NFW" # Primo ingrediente: profilo di densità dell'alone di materia oscura
# Navarro - Frenk - White, lo standard più comune
spettro_noBAO.bias_model = 'ST99' # Secondo ingrediente, il "bias dell'alone"
spettro_noBAO.hmf_model = 'ST' # Terzo ingrediente: la funzione di massa, quanti (in numero) aloni abbiamo nell'universo per ogni valore di massa dell'alone?
Modello basato su parametri cosmologici alternativi, in particolare, diverso $H_{0}$, diverso $\Omega_{\rm{materia \ oscura},0}$
pk_contrasto_BAO = DMHaloModel()
pk_contrasto_BAO.z = hm_ref.z
pk_contrasto_BAO.tracer_profile_model = "NFW" # Primo ingrediente: profilo di densità dell'alone di materia oscura
# Navarro - Frenk - White, lo standard più comune
pk_contrasto_BAO.bias_model = 'ST99' # Secondo ingrediente, il "bias dell'alone"
pk_contrasto_BAO.hmf_model = 'ST' # Terzo ingrediente: la funzione di massa, quanti (in numero) aloni abbiamo nell'universo per ogni valore di massa dell'alone?
pk_contrasto_BAO.cosmo_params = {
'H0': 67.0,
'Om0':0.15
}
pk_contrasto_BAO.sigma_8 = 0.68#spettro_BAO.sigma_8#0.68
pk_contrasto_noBAO = DMHaloModel()
redshift = 0.
pk_contrasto_noBAO.z = redshift
pk_contrasto_noBAO.tracer_profile_model = "NFW" # Primo ingrediente: profilo di densità dell'alone di materia oscura
# Navarro - Frenk - White, lo standard più comune
pk_contrasto_noBAO.bias_model = 'ST99' # Secondo ingrediente, il "bias dell'alone"
pk_contrasto_noBAO.hmf_model = 'ST' # Terzo ingrediente: la funzione di massa, quanti (in numero) aloni abbiamo nell'universo per ogni valore di massa dell'alone?
pk_contrasto_noBAO.cosmo_params = {
'H0': 75.0,
'Om0':0.15,
}
pk_contrasto_noBAO.transfer_model = 'EH_NoBAO'
pk_contrasto_noBAO.sigma_8 = 0.68#spettro_BAO.sigma_8#0.68
Fin qui, abbiamo generato i dati, le nostre "osservazioni" del mondo fisico. Questi però sono viziati da incertezze, che possono dipendere da effetti strumentali o astrofisici sottraibili a posteriori, noti e solitamente sottratti o quantificati a parte, gli errori sistematici, e residui errori di misura che, nel caso della cosmologia, combinano proprietà statistiche dei dati e delle misure e proprietà fisiche (con un pizzico in più di statistica nel caso della cosmologia, ricordiamo che lo spettro di potenza è una misura di una correlazione statistica!). Per semplicità, le incertezze che creo nella prossima cella non sono quelle che effettivamente si assegnano alle osservazioni, ma costruite manualmente e assolutamente "finte", solo per poterci poi ragionare sopra.
Ciò che è fondamentale tenere a mente è però questo concetto: sono le dimensioni delle incertezze che ci consentono di distinguere ("risolvere") tra diversi modelli fisici, perché più queste sono grandi più, intuitivamente, diversi modelli saranno compatibili ("fitteranno") i dati. (Nota: to fit, in inglese, vuol dire appunto essere adatto, appropriato, in forma!)
# Genero delle finte incertezze sui dati, per riprodurre il tipico andamento a fisarmonica ho trovato questa combinazione
uncert = hm_ref.power_auto_matter_fnc(kbins)**0.7*kbins + hm_ref.power_auto_matter_fnc(kbins)*np.exp(-kbins*30)
plt.figure(figsize=(18,12))
plt.errorbar(kbins, p_k,yerr= uncert,label='dati',ls='')
plt.loglog(hm_ref.k_hm,spettro_BAO.power_auto_matter_fnc(hm_ref.k_hm),label='cosmologia di riferimento con BAO')
plt.loglog(hm_ref.k_hm,spettro_noBAO.power_auto_matter_fnc(hm_ref.k_hm),label='cosmologia di riferimento no BAO')
plt.loglog(hm_ref.k_hm,pk_contrasto_BAO.power_auto_matter_fnc(hm_ref.k_hm),label='cosmologia alternativa con BAO')
plt.loglog(hm_ref.k_hm,pk_contrasto_noBAO.power_auto_matter_fnc(hm_ref.k_hm),label='cosmologia alternativa no BAO')
plt.xscale('log')
plt.yscale('log')
plt.xlabel("k [$Mpc^{-1} h$]",size=32)
plt.ylabel(r" P(k)$ \ [{ Mpc^3}h^{-3}]$",size=32);
plt.legend()
plt.grid()
#plt.xlim(0.008,1)
/home/andrea/anaconda3/envs/progetto_1/lib/python3.11/site-packages/halomod/tools.py:653: UserWarning: You are using an un-normalized mass function and bias function pair.Bias ST99 has the following paired HMF model: (<class 'hmf.mass_function.fitting_functions.SMT'>,). Matter correlations are not well-defined. warnings.warn(
Che cosa si può dire dei vari modelli? Sono entrambi affidabili?
dati = p_k[cut_index_in:cut_index_fin+1]
incertezze = uncert[cut_index_in:cut_index_fin+1]
plt.figure(figsize=(10,5))
# Vediamo il nostro rapporto segnale rumore
plt.title('Rapporto segnale rumore (SNR)')
plt.semilogx(kbins[cut_index_in:cut_index_fin+1],dati/incertezze)
plt.xlabel("k [$Mpc^{-1} h$]",size=32)
plt.grid()
Qui comincia il nostro fit dei dati: abbiamo le misure, abbiamo due modelli di spettro di potenza, entrambi costruiti in questo modo: assumiamo che ci sia una componente "garantita" dello spettro di potenza senza oscillazioni, $P_{no BAO}$ e ad essa aggiungiamo un termine con le oscillazioni, in questo modo: $$P_\text{modello} (k) = A ( P_{BAO}(k/\alpha) - P_{no BAO}(k/\alpha) ) + P_{no BAO(k)} + B * k $$
Allora:
Il fit si fa in questo modo:
#Fit con due modelli
def lnlike_1(theta):
# Modello "vero"
A, alpha , B = theta
modello = np.zeros(len(k_data))
funzione_fit = lambda k: A * ( spettro_BAO.power_auto_matter_fnc(k/alpha) - spettro_noBAO.power_auto_matter_fnc(k/alpha) ) + spettro_noBAO.power_auto_matter_fnc(k) + B * k
modello = funzione_fit(k_data)
scarti_su_incertezze = np.sum(( ( dati - modello) / incertezze ) ** 2)
return -0.5 * scarti_su_incertezze
def lnlike_2(theta):
# Modello "di contrasto"
A, alpha ,B = theta
modello = np.zeros(len(k_data)) #monopole model
funzione_fit = lambda k: A * ( pk_contrasto_BAO.power_auto_matter_fnc(k) - pk_contrasto_noBAO.power_auto_matter_fnc(k/alpha) ) + pk_contrasto_noBAO.power_auto_matter_fnc(k) + B * k
modello = funzione_fit(k_data)
scarti_su_incertezze = np.sum(( ( dati - modello) / incertezze ) ** 2)
return -0.5*scarti_su_incertezze
"Fissiamo i prior".
Esistono due modi fondamentali per fare analisi statistiche:
In generale, tecniche di combinazione tra le osservabili, scelta dei prior in modo che siano vincolanti, ad esempio ricorrendo alle osservazioni di Planck della radiazione cosmica di fondo, rivestono un ruolo importantissimo nell'analisi dei dati in fisica.
def lnprob_1(theta):
lp = lnprior_1(theta)
if not np.isfinite(lp):
return -np.inf
return lp + lnlike_1(theta)
def lnprior_1(theta):
A ,alpha, B = theta
if 0 < A < 2 and 0. < alpha < 2 and 0 < B <1:
return 0.0
return -np.inf
def lnprob_2(theta):
lp = lnprior_2(theta)
if not np.isfinite(lp):
return -np.inf
return lp + lnlike_2(theta)
def lnprior_2(theta):
A, alpha ,B = theta
if 0 < A < 2 and 0. < alpha < 2 and 0 < B <1:
return 0.0
return -np.inf
# Eseguiamo i nostri due fit
import emcee # importiamo il programma che farà il lavoro per noi
nwalkers = 100 # imponiamo un certo numero di punti di partenza all'esplorazione dello spazio dei parametri
niter=10000 # imponiamo un certo numero di passi per assicurarci che la catena giunga poi effettivamente a convergenza
# e non "si inganni" in un minimo locale provvisorio.
npars = 3 #numero di parametri
fiduc=np.array([1,1,0]) #valori fiduciari, quelli che ci aspettiamo saranno il valore di "Best fit"
ndim = len(fiduc)
randy = fiduc + 1e-2 * np.random.randn(nwalkers,ndim) # comincia ad esplorare ad una certa distanza dal valore di partenza
pos = fiduc + randy
from multiprocessing import Pool # per fare prima, parallelizziamo l'operazione, distribuendola sui vari processori del
# computer
with Pool() as pool:
sampler_1=emcee.EnsembleSampler(nwalkers,ndim,lnprob_1,pool=pool)
sampler_1.run_mcmc(pos, niter, progress=True);
0%| | 0/10000 [00:00<?, ?it/s]/home/andrea/anaconda3/envs/progetto_1/lib/python3.11/site-packages/emcee/moves/red_blue.py:99: RuntimeWarning: invalid value encountered in scalar subtract lnpdiff = f + nlp - state.log_prob[j] 100%|████████████████████████████████████████████████████████████████████| 10000/10000 [01:23<00:00, 120.27it/s]
from multiprocessing import Pool
with Pool() as pool:
sampler_2=emcee.EnsembleSampler(nwalkers,ndim,lnprob_2,pool=pool)
sampler_2.run_mcmc(pos, niter, progress=True);
100%|████████████████████████████████████████████████████████████████████| 10000/10000 [01:21<00:00, 123.09it/s]
La Monte Carlo Markov Chain (MCMC) non fornisce un valore di best fit per ogni parametro (come per l'approccio inferenziale), ma una distribuzione di probabilità a posteriori (il posterior) che, di solito, è gaussiana, con un massimo ben definito che rappresenta il valore di best-fit. Ogni valore, in teoria, del posterior, è ammesso, ma il massimo è favorito in quanto più probabile. D'altra parte, non è detto che i posterior siano così gaussiani come ci aspetteremmo...
nburn = 400
labels = [r'A',r'$\alpha$','B']
flat_samples_1 = sampler_1.get_chain(discard=int(nburn*2), thin=int(0.5*nburn), flat=True)
names = ['A','alpha','B']
markers={'A':fiduc[0],'$\alpha$':fiduc[1],'B':fiduc[2],}
samples_1 = MCSamples(samples=flat_samples_1,names = names, labels = labels)
g = plots.get_subplot_plotter()
g.settings.figure_legend_frame = False
g.settings.alpha_filled_add=0.4
g.settings.title_limit_fontsize = 14
g.triangle_plot([samples_1], filled=True,title_limit=None,markers=markers)
WARNING:root:auto bandwidth for alpha very small or failed (h=0.0014493129346851159,N_eff=4600.0). Using fallback (h=0.0017232741113171522)
Removed no burn in
nburn = 400
labels = [r'A',r'\alpha','B']
flat_samples_2 = sampler_2.get_chain(discard=int(nburn*2), thin=int(0.5*nburn), flat=True)
names = ['A','alpha','B']#'npar_fid','nperp_fid']#,'k_par_quadru','k_perp_quadru']
markers={'A':fiduc[0],'$\alpha$':fiduc[1],'B':fiduc[2],}#'npar_fid':npar_fid,'nperp_fid':nperp_fid}#,'k_par_quadru':k_par_fid_quadru/kpar0,'k_perp_quadru':k_perp_fid_quadru/kperp0}
#markers={bg_fid,'bh':pref_H_fid,'r':R_fid,'sigma':sigmar_fid,'k_par_mono':k_par_fid_mono/kpar0,'k_perp_mono':k_perp_fid_mono/kperp0,'k_par_quadru':k_par_fid_quadru/kpar0,'k_perp_quadru':k_perp_fid_quadru/kperp0}
samples_2 = MCSamples(samples=flat_samples_2,names = names, labels = labels)
g = plots.get_subplot_plotter()
g.settings.figure_legend_frame = False
g.settings.alpha_filled_add=0.4
g.settings.title_limit_fontsize = 14
g.triangle_plot([samples_2], filled=True,title_limit=None,markers=markers)
WARNING:root:auto bandwidth for A very small or failed (h=0.001748086484650311,N_eff=4600.0). Using fallback (h=0.0020245372894739525) WARNING:root:auto bandwidth for alpha very small or failed (h=0.0011030819747355634,N_eff=4600.0). Using fallback (h=0.00206712766195534)
Removed no burn in
Siamo principalmente interessati ad $\alpha$: confrontiamo tra i due risultati e proviamo a commentarli
# 1D marginalized comparison plot
g = plots.get_single_plotter(width_inch=5)
g.plot_1d([samples_1, samples_2], 'alpha',texstize=30)
WARNING:root:auto bandwidth for alpha very small or failed (h=0.0014493129346851159,N_eff=4600.0). Using fallback (h=0.0017232741113171522) WARNING:root:auto bandwidth for alpha very small or failed (h=0.0011030819747355634,N_eff=4600.0). Using fallback (h=0.00206712766195534)
Che cosa ci dicono questi parametri? Sono compatibili? Sono in tensione?
from IPython.display import display, Math
sample1=np.zeros(len(fiduc))
uncert_par_1=np.zeros([ndim,2]) #store +- uncertainties
for i in range(ndim):
mcmc = np.percentile(flat_samples_1[:, i], [16, 50, 84]) #best-fit value =50th percentile
sample1[i]=mcmc[1]#round(mcmc[1], 2) samples.means[i]
q = np.diff(mcmc)
uncert_par_1[i]=q
q = np.diff(mcmc)
txt = "\mathrm{{{3}}} = {0:.3f}_{{-{1:.3f}}}^{{{2:.3f}}}"
txt = txt.format(mcmc[1], q[0], q[1], labels[i])
display(Math(txt))
sample2=np.zeros(len(fiduc))
uncert_par_2=np.zeros([ndim,2]) #store +- uncertainties
for i in range(ndim):
mcmc = np.percentile(flat_samples_2[:, i], [16, 50, 84]) #best-fit value =50th percentile
sample2[i]=mcmc[1]#round(mcmc[1], 2) samples.means[i]
q = np.diff(mcmc)
uncert_par_2[i]=q
q = np.diff(mcmc)
print(q)
txt = "\mathrm{{{3}}} = {0:.3f}_{{-{1:.3f}}}^{{{2:.3f}}}"
txt = txt.format(mcmc[1], q[0], q[1], labels[i])
display(Math(txt))
[0.01 0.01]
[0.007 0.008]
[0.367 0.258]
$\alpha$, specialmente il secondo, non è molto gaussiano, le due incertezze non sono molto simmetriche. Assumiamo che valga l'incertezza maggiore, e calcoliamo allora questa quantità $\dfrac{ \alpha_1 - \alpha_2 }{ \sqrt{ \sigma_{1}^2 + \sigma_{2}^2 }} $
Che ci dice quanto, e se, i due parametri sono in tensione.
(sample1[1]-sample2[1])/np.sqrt(np.min(uncert_par_1[1])**2 + np.min(uncert_par_2[1])**2) # rivedere, ho un po' "adattato"
53.69724622820244
I due parametri sono in tensione a 53 $\sigma$, decisamente incompatibili! Ma a sua volta, il risultato sarà compatibile con la cosmologia fiduciaria, potremo dire che $\alpha$ vale 1 entro le incertezze?
(sample1[1]-1)/np.sqrt(np.max(uncert_par_1[1])**2)
0.9697161459348516
Sì, perché la differenza è di circa 1 sigma (non un fit perfetto, sicuramente) ma minore dei tipici criteri di 3 o 5 sigma adottati in fisica per stabilire la significatività di una differenza, o di una compatibilità, tra i dati e il modello.
# Confrontiamo i due risultati
funzione_fit_1 = lambda k: sample1[0] * ( spettro_BAO.power_auto_matter_fnc(k/sample1[1]) - spettro_noBAO.power_auto_matter_fnc(k/sample1[1]) ) + spettro_noBAO.power_auto_matter_fnc(k) + sample1[2] * k
funzione_fit_2 = lambda k: sample2[0] * ( pk_contrasto_BAO.power_auto_matter_fnc(k/sample2[1]) - pk_contrasto_noBAO.power_auto_matter_fnc(k/sample2[1]) ) + pk_contrasto_noBAO.power_auto_matter_fnc(k) + sample2[2] * k
plt.figure(figsize=(10,6))
plt.errorbar(k_data, dati,yerr= incertezze,label='dati')
#plt.loglog(hm.k_hm,P_nwh(hm.k_hm))
#plt.loglog(hm.k_hm,P_wh(hm.k_hm))
plt.loglog(k_data,funzione_fit_1(k_data),label='modello 1')
plt.loglog(k_data,funzione_fit_2(k_data),label='modello 2')
plt.xscale('log')
plt.yscale('log')
plt.legend()
plt.xlabel("k [$Mpc^{-1} h$]")
plt.ylabel(r"$\rm P(k) \ [{\rm Mpc^3}h^{-3}]$");
plt.xlim(2*1e-2,1)
plt.grid()
Peter Collier, A most incomprehensible thing: introduzione alla relatività generale, autopubblicata su Amazon, comprensiva di cosmologia. Molto interessante, tra il qualitativo ed il qualitativo, introduzione ben fatta per appassionati che, come nel caso del libro di Penrose, desiderino impadronirsi delle basi del linguaggio matematico della relatività generale, "leggere le formule" e capire come sono fatte, che cosa significhino e quali siano le incognite delle equazioni che i fisici vogliono misurare.
Edward W. Kolb, Michael S. Turner, The Early Universe, anche in questo caso disponibile su Amazon come ristampa, solo in inglese. Dall'universo omogeneo alle basi della cosmologia quantistica, un testo fondamentale specialmente per capire l'evoluzione dell'universo nella sua fase primordiale e "calda", quando vengono sintetizzati gli elementi costitutivi della materia che conosciamo (bariogenesi e nucleosintesi primordiale), che copre però anche molto bene la teoria delle perturbazioni e la formazione delle strutture.
Attilio Ferrari, Stelle, Galassie e Universo. Fondamenti di Astrofisica, Springer Verlag, 2011. Manuale universitario di astrofisica e cosmologia (con utile appendice online anche da scaricare), redatto da uno dei più importanti docenti di astrofisica dell'Università di Torino, radioastronomo di formazione, attento alla divulgazione e a lungo presidente del Planetario e Museo dello spazio di Torino, nonché direttore dell'osservatorio astronomico di Pino Torinese.
John Peacock, Cosmological Physics, Cambridge University Press, 1998. Manuale universitario, copre in realtà molto di più della cosmologia (ma, d'altra parte, la cosmologia è per definizione una branca in cui tutta la conoscenza accumulata dalla fisica nel tempo diventa in qualche modo necessaria) di quanto il titolo lasci intuire, dalla relatività generale alla teoria dei campi all'astrofisica. Ha un taglio un po' datato, più astrofisico e meno statistico, ma proprio per questo può aiutare a cogliere il senso "osservativo" di certi concetti.
Scott Dodelson, Fabian Schmidt, Modern Cosmology, 2020. Semplicemente IL MANUALE di cosmologia, per impararne davvero il linguaggio contemporaneo, basato sulla teoria perturbativa e l'inferenza bayesiana. Esiste anche una prima edizione, ancora disponibile pubblicata nel 2003. Un poco datata, non altrettanto onnicomprensiva, ma validissima allo stesso modo. Si tratta di un saggio universitario gradevole, ma richiede già una preparazione avanzata, è oggetivamente difficile per un liceale curioso. Per chi desiderasse capire "come funziona" la cosmologia, a cavallo tra quantitativo e qualitativo, i testi precedenti restano quelli consigliati.
Note di Alessandro Marconi per il corso di Cosmologia all'Università di Firenze. Si tratta di un valido manuale di cosmologia in italiano
che copre tutto quanto è necessario sapere per avere una formazione sull'argomento, dall'espansione dell'universo alla CMB. Ha un taglio più vicino al Peacock che al Dodelson, lievemente più statistico che astrofisico, ma sicuramente tra i migliori testi di riferimento a disposizione in italiano. Ci sono spiegazioni "in prosa", grafici ed altri materiali utili anche per un liceale, se queste lezioni hanno suscitato passione ed interesse per l'argomento, magari un giorno tornerà utile.